﻿#include <iostream>
#include <iomanip>

#include <Eigen/Sparse>
#include <Eigen/SparseLU>
#include <fstream>
#include <time.h>

#include "gmsh.h"
#include "point.h"
#include "Triangle.h"

using namespace std;

//example index
int ex_index = 3;

typedef Eigen::SparseMatrix<double> SpMat;//declares a sparse matrix type of double
typedef Eigen::VectorXd Rsv;//declares a right side value type
typedef Eigen::Triplet<double> Tri;//declares a triplet to stores the matrix index and values, then convert it to sparse matrix(with high performance)

double analytical_solution(Point p, int der_x, int der_y)
{
	//example 1 in chapter 3
	double x = p.getX();
	double y = p.getY();

	if (ex_index == 1)
	{
		if (der_x == 0 && der_y == 0)
		{
			return (p.getX()*p.getY()*(1 - p.getX() / 2.0)*(1 - p.getY())*exp(p.getX() + p.getY()));

		}
		else if (der_x == 1 && der_y == 0)
		{
			return x * y*(1 - x / 2)*(1 - y)*exp(x + y) - x * y*(1 - y)*exp(x + y) / 2 + y * (1 - x / 2)*(1 - y)*exp(x + y);
		}
		else if (der_x == 0 && der_y == 1)
		{
			return x * y*(1 - x / 2)*(1 - y)*exp(x + y) - x * y*(1 - x / 2)*exp(x + y) + x * (1 - x / 2)*(1 - y)*exp(x + y);
		}
		else if (der_x == 1 && der_y == 1)
		{
			return (x*y*(x - 2)*(y - 1) + x * y*(x - 2) + x * y*(y - 1) + x * y + x * (x - 2)*(y - 1) + x * (y - 1) + y * (x - 2)*(y - 1) + y * (x - 2) + (x - 2)*(y - 1))*exp(x + y) / 2;
		}
	}

	else if (ex_index == 2 || ex_index == 3)
	{
		//example 2 and 3 in chapter 3
		if (der_x == 0 && der_y == 0)
		{
			return exp(x + y);

		}
		else if (der_x == 1 && der_y == 0)
		{
			return exp(x + y);
		}
		else if (der_x == 0 && der_y == 1)
		{
			return exp(x + y);
		}
		else if (der_x == 1 && der_y == 1)
		{
			return exp(x + y);
		}
		else
			return 0;
	}
	return 0;


}
//coefficients
double coe_function(Point p)
{
	return 1.0;
}
double load_function(Point p)
{
	if (1 == ex_index)
	{
		//example 1 in chapter 3
		return  -(p.getY()*(1 - p.getY())*(1 - p.getX() - pow(p.getX(), 2) / 2.0)*exp(p.getX() + p.getY())
			+ p.getX()*(1 - p.getX() / 2.0)*(-3 * p.getY() - pow(p.getY(), 2))*exp(p.getX() + p.getY()));
	}
	else if (2 == ex_index || 3 == ex_index)
	{
		//example 2 in chapter 3
		return -2.0*exp(p.getX() + p.getY());
	}
	return 0;

}
double left_bound_value(Point p)
{
	if (1 == ex_index)
	{
		//example 1 in chapter 3
		return -1.5*p.getY()*(1 - p.getY())*exp(-1.0 + p.getY());
	}

	else if (2 == ex_index || 3 == ex_index)
	{
		//example 2 in chapter 3
		return exp(-1 + p.getY());
	}
	return 0;

}
double right_bound_value(Point p)
{
	if (1 == ex_index)
	{
		//example 1 in chapter 3
		return 0.5*p.getY()*(1 - p.getY())*exp(1 + p.getY());
	}

	else if (2 == ex_index || 3 == ex_index)
	{
		//example 2 in chapter 3
		return exp(1 + p.getY());
	}
	return 0;

}
double bottom_bound_value(Point p)
{
	if (1 == ex_index)
	{
		//example 1 in chapter 3
		return -2 * p.getX()*(1 - p.getX() / 2.0)*exp(p.getX() - 1);
	}

	else if (2 == ex_index )
	{
		//example 2 in chapter 3
		return -exp(p.getX() - 1);
	}
	else if (3==ex_index)
	{
		return 0;
	}
	return 0;

}
double up_bound_value(Point p)
{
	if (1 == ex_index)
	{
		//example 1 in chapter 3
		return 0;
	}

	else if (2 == ex_index || 3 == ex_index)
	{
		//example 2 in chapter 3
		return exp(p.getX() + 1);
	}
	return 0;

}

//nalba_u+ru=q
double r(Point x)
{
	return 1;
}

//The mesh generated by gambit or gmsh is Pb and Tb
//so we need get P and T from Pb and Tb
bool generate_P_T(BASIS_TYPE basis_type_trial, BASIS_TYPE basis_type_test, vector<Point> &Nodes, vector<Triangle> &Triangles,
	vector<Point> &Pb_trial, vector<Triangle> &Tb_trial,
	vector<Point> Pb_test, vector<Triangle> Tb_test)
{

	//P and T
	if (basis_type_test == TWO_DIM_LINEAR)//basis function type
	{
		//For linear basis function, P and T are Pb and Tb
		Nodes = Pb_test;
		Triangles = Tb_test;

	}
	else if (basis_type_test == TWO_DIM_BILINEAR)
	{
		Nodes = Pb_test;
		Triangles = Tb_test;
	}
	else if (basis_type_test == TWO_DIM_QUADRATIC)
	{
		for (unsigned int i = 0; i < Tb_test.size(); ++i)
		{
			//The P has duplicate nodes stored, but it doesn't affect the computation because we only need the right coordinates
			Triangle T;
			for (int j = 0; j < Tb_test[0].get_ele_node_num(); j = j + 2)
			{
				Nodes.push_back(Pb_test[Tb_test[i].getIds()[j]]);
				T.append(j / 2 + i * 3);
			}
			Triangles.push_back(T);
		}
	}
	else if (basis_type_test == TWO_DIM_BIQUADRATIC)
	{
		for (unsigned int i = 0; i < Tb_test.size(); ++i)
		{
			//The P has duplicate nodes stored, but it doesn't affect the computation because we only need the right coordinates
			Triangle T;
			for (int j = 0; j < Tb_test[0].get_ele_node_num() - 1; j = j + 2)
			{
				Nodes.push_back(Pb_test[Tb_test[i].getIds()[j]]);
				T.append(j / 2 + i * 4);
			}
			Triangles.push_back(T);
		}
	}

	//Pb_trial and Tb_trial
	Pb_trial = Pb_test;
	Tb_trial = Tb_test;

	return true;
}

bool export_to_tecplot(vector<Point> Nodes, vector<Triangle> Triangles, string filename, Rsv x, BASIS_TYPE basis_type)
{
	ofstream f(filename);
	f << "variables=x,y,z,s_fem,s_ana" << endl;
	if (TWO_DIM_LINEAR == basis_type || TWO_DIM_QUADRATIC == basis_type)
	{
		f << "ZONE N=\t" << Nodes.size() << ",\t" << "E=\t" << Triangles.size() << ",\tDATAPACKING=POINT" << ",\tZONETYPE=\t FETRIANGLE" << endl;
	}
	else if (TWO_DIM_BILINEAR == basis_type || TWO_DIM_BIQUADRATIC == basis_type)
	{
		f << "ZONE N=\t" << Nodes.size() << ",\t" << "E=\t" << Triangles.size() << ",\tDATAPACKING=POINT" << ",\tZONETYPE=\t FEQUADRILATERAL" << endl;
	}
	for (unsigned int i = 0; i < Nodes.size(); ++i)
	{
		f << Nodes[i].getX() << "\t" << Nodes[i].getY() << "\t" << Nodes[i].getZ() << "\t" << x(i) << "\t" << analytical_solution(Nodes[i], 0, 0) << endl;
	}
	for (unsigned int i = 0; i < Triangles.size(); ++i)
	{
		if (TWO_DIM_LINEAR == basis_type)
		{
			for (unsigned int j = 0; j < 3; ++j)
			{
				f << Triangles[i].getIds()[j] + 1 << "\t";
			}
		}
		else if (TWO_DIM_QUADRATIC == basis_type)
		{
			for (unsigned int j = 0; j < 3; ++j)
			{
				f << Triangles[i].getIds()[2 * j] + 1 << "\t";
			}
		}
		else if (TWO_DIM_BILINEAR == basis_type)
		{
			for (unsigned int j = 0; j < 4; ++j)
			{
				f << Triangles[i].getIds()[j] + 1 << "\t";
			}
		}
		else if (TWO_DIM_BIQUADRATIC == basis_type)
		{
			for (unsigned int j = 0; j < 4; ++j)
			{
				f << Triangles[i].getIds()[2 * j] + 1 << "\t";
			}

		}

		f << endl;

	}
	f.close();
	return true;

}
//Tecplot does not support the quadratic element, so we export results to vtk and visualizes them with Paraview 
bool export_to_vtk(vector<Point>Pb, vector<Triangle> Tb, string filename, Rsv x)
{
	map<ELE_TYPE, int> type_ele_vtk = {
		{ELE_3,5},
		{ELE_4,9},
		{ELE_6,22},
		{ELE_9,28}
	};
	ofstream f(filename);
	f << "# vtk DataFile Version 3.0" << endl;
	f << "Result" << endl;
	f << "ASCII" << endl;
	f << "DATASET UNSTRUCTURED_GRID" << endl;

	f << "POINTS \t" << Pb.size() << "\tdouble" << endl;
	for (unsigned int i = 0; i < Pb.size(); ++i)
	{
		f << Pb[i].getX() << "\t" << Pb[i].getY() << "\t" << Pb[i].getZ() << endl;
	}
	f << "CELLS \t " << Tb.size() << "\t" << (Tb[0].getIds().size() + 1)* Tb.size() << endl;
	for (unsigned int i = 0; i < Tb.size(); ++i)
	{

		f << Tb[i].getIds().size() << "\t";

		if (Tb[i].get_ele_type() == ELE_3 || Tb[i].get_ele_type() == ELE_4)
		{
			for (unsigned int j = 0; j < Tb[i].getIds().size(); ++j)
			{
				f << Tb[i].getIds()[j] << "\t";
			}
		}
		else if (Tb[i].get_ele_type() == ELE_6)
		{
			for (unsigned int j = 0; j < Tb[i].getIds().size(); j += 2)
			{
				f << Tb[i].getIds()[j] << "\t";
			}
			for (unsigned int j = 1; j < Tb[i].getIds().size(); j += 2)
			{
				f << Tb[i].getIds()[j] << "\t";
			}
		}
		else if (Tb[i].get_ele_type() == ELE_9)
		{
			for (unsigned int j = 0; j < Tb[i].getIds().size() - 1; j += 2)
			{
				f << Tb[i].getIds()[j] << "\t";
			}
			for (unsigned int j = 1; j < Tb[i].getIds().size() - 1; j += 2)
			{
				f << Tb[i].getIds()[j] << "\t";
			}
			f << Tb[i].getIds()[8] << "\t";
		}
		f << endl;

	}
	f << "CELL_TYPES\t" << Tb.size() << endl;

	for (unsigned int i = 0; i < Tb.size(); ++i)
	{
		f << type_ele_vtk[Tb[i].get_ele_type()] << endl;
	}

	f << "POINT_DATA\t" << Pb.size() << endl;

	f << "SCALARS s_fem double 1" << endl;

	f << "LOOKUP_TABLE default" << endl;
	for (unsigned int i = 0; i < Pb.size(); ++i)
	{

		f << x(i) << endl;

	}
	f << "SCALARS s_ana double 1" << endl;
	f << "LOOKUP_TABLE default" << endl;

	for (unsigned int i = 0; i < Pb.size(); ++i)
	{

		f << analytical_solution(Pb[i], 0, 0) << endl;

	}
	f.close();
	return true;


}

//Generate mesh with gmsh
//input:
//left down corner point and right up corner point
//return:
//all mesh nodes (Points) and elements(Triangles) and all the boundary nodes and boundary edges
bool generate_mesh_with_gmsh(double x_left_down, double y_left_down, double x_right_up, double y_right_up, double mesh_size,
	vector<Point> &Nodes, vector<Triangle> &triangles, vector<vector<boundary_node>> &allBoundNodes,
	vector<vector<boundary_edge>> &allBoundEdges, ELE_TYPE &ele_type)
{
	ele_type = ELE_3;
	if (x_right_up <= x_left_down || y_right_up <= y_left_down)
	{
		cout << "\nerror:Wrong input for geometry, can't generate mesh\n";
		return false;
	}

	gmsh::initialize();

	gmsh::model::add("test");

	//Add point
	int p1 = gmsh::model::geo::addPoint(x_left_down, y_left_down, 0);
	int p2 = gmsh::model::geo::addPoint(x_right_up, y_left_down, 0);
	int p3 = gmsh::model::geo::addPoint(x_right_up, y_right_up, 0);
	int p4 = gmsh::model::geo::addPoint(x_left_down, y_right_up, 0);

	//Add line
	int l1 = gmsh::model::geo::addLine(p1, p2);
	int l2 = gmsh::model::geo::addLine(p2, p3);
	int l3 = gmsh::model::geo::addLine(p3, p4);
	int l4 = gmsh::model::geo::addLine(p4, p1);

	//Add physical group for boundary
	int dim = 1;
	int boundaryTagBottom = gmsh::model::addPhysicalGroup(dim, { l1 });
	gmsh::model::setPhysicalName(dim, boundaryTagBottom, "Bottom");

	int boundaryTagRight = gmsh::model::addPhysicalGroup(dim, { l2 });
	gmsh::model::setPhysicalName(dim, boundaryTagBottom, "Right");

	int boundaryTagUp = gmsh::model::addPhysicalGroup(dim, { l3 });
	gmsh::model::setPhysicalName(dim, boundaryTagBottom, "Up");

	int boundaryTagLeft = gmsh::model::addPhysicalGroup(dim, { l4 });
	gmsh::model::setPhysicalName(dim, boundaryTagBottom, "Left");

	//Add curve loop
	int c1 = gmsh::model::geo::addCurveLoop({ l1,l2,l3,l4 });

	//Add surface
	int s1 = gmsh::model::geo::addPlaneSurface({ c1 });

	dim = 2;
	gmsh::model::addPhysicalGroup(dim, { c1 });

	//synchronize
	gmsh::model::geo::synchronize();

	//define the mesh size
	gmsh::model::mesh::field::add("Box", 1);
	gmsh::model::mesh::field::setNumber(1, "VIn", mesh_size);
	gmsh::model::mesh::field::setNumber(1, "VOut", mesh_size);
	gmsh::model::mesh::field::setNumber(1, "XMin", x_left_down);
	gmsh::model::mesh::field::setNumber(1, "XMax", x_right_up);
	gmsh::model::mesh::field::setNumber(1, "YMin", y_left_down);
	gmsh::model::mesh::field::setNumber(1, "YMax", y_right_up);
	gmsh::model::mesh::field::setAsBackgroundMesh(1);
	gmsh::option::setNumber("Mesh.MeshSizeExtendFromBoundary", 0);
	gmsh::option::setNumber("Mesh.MeshSizeExtendFromBoundary", 0);

	gmsh::model::mesh::generate(2);

	//gmsh::write("demo.msh");

	//gmsh::fltk::run();

	//Get nodes
	dim = -1;
	int tag = -1;
	vector<double> coord, parametricCoord;
	vector<size_t> nodeTags;
	gmsh::model::mesh::getNodes(nodeTags, coord, parametricCoord, dim, tag);
	cout << "The number of nodes: " << coord.size() / 3 << endl;

	//nodes.reserve(1000);
	Point p;
	for (unsigned int i = 0; i < coord.size(); i += 3)
	{
		double x = coord[i];
		double y = coord[i + 1];
		double z = coord[i + 2];
		p = Point(x, y, z);
		Nodes.push_back(p);
	}

	//Get elements	
	vector<int> elementTypes;
	vector<vector<size_t> > elementTags, nodeTagss;
	gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTagss);
	cout << "The number of elements: " << nodeTagss[1].size() / 3 << endl;

	unsigned int id1, id2, id3;

	for (unsigned int i = 0; i < nodeTagss[1].size(); i += 3)
	{
		Triangle tri;
		id1 = nodeTagss[1][i] - 1;
		id2 = nodeTagss[1][i + 1] - 1;
		id3 = nodeTagss[1][i + 2] - 1;
		tri.append(id1);
		tri.append(id2);
		tri.append(id3);
		triangles.push_back(tri);
	}

	cout << "get boundary nodes" << endl;
	//get boundary nodes
	dim = 1;
	vector<size_t> bouNodesIdsLeft, bouNodesIdSRight, bouNodesIdsBottom, BouNodesIdsUp;
	gmsh::model::mesh::getNodesForPhysicalGroup(dim, boundaryTagLeft, bouNodesIdsLeft, coord);
	gmsh::model::mesh::getNodesForPhysicalGroup(dim, boundaryTagBottom, bouNodesIdsBottom, coord);
	gmsh::model::mesh::getNodesForPhysicalGroup(dim, boundaryTagRight, bouNodesIdSRight, coord);
	gmsh::model::mesh::getNodesForPhysicalGroup(dim, boundaryTagUp, BouNodesIdsUp, coord);


	vector<boundary_node> boundaryNodesIdsLeft(bouNodesIdsLeft.size());
	// just because index from start 0
	for (unsigned int i = 0; i < bouNodesIdsLeft.size(); i++)
	{
		boundaryNodesIdsLeft[i].glo_node_index = (bouNodesIdsLeft[i] - 1);
	}
	allBoundNodes.push_back(boundaryNodesIdsLeft);

	vector<boundary_node> boundaryNodesIdsBottom(bouNodesIdsBottom.size());
	// just because index from start 0
	for (unsigned int i = 0; i < bouNodesIdsBottom.size(); i++)
	{
		boundaryNodesIdsBottom[i].glo_node_index = (bouNodesIdsBottom[i] - 1);
	}
	allBoundNodes.push_back(boundaryNodesIdsBottom);


	vector<boundary_node> boundaryNodesIdsRight(bouNodesIdSRight.size());
	// just because index from start 0
	for (unsigned int i = 0; i < bouNodesIdSRight.size(); i++)
	{
		boundaryNodesIdsRight[i].glo_node_index = (bouNodesIdSRight[i] - 1);
	}
	allBoundNodes.push_back(boundaryNodesIdsRight);


	vector<boundary_node> boundaryNodesIdsUp(BouNodesIdsUp.size());
	// just because index from start 0
	for (unsigned int i = 0; i < BouNodesIdsUp.size(); i++)
	{
		boundaryNodesIdsUp[i].glo_node_index = (BouNodesIdsUp[i] - 1);
	}
	allBoundNodes.push_back(boundaryNodesIdsUp);


	//get boundary edges
	//bottom boundary
	cout << "get boundary edges" << endl;
	dim = 1;
	int dim2 = 2;
	bool strict = true;
	tag = l1;
	vector<size_t> elementTag;
	vector<boundary_edge> bottom_bound_egdes;
	gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTagss, dim, tag);
	for (unsigned int k = 0; k < nodeTagss[0].size(); k += 2)
	{
		id1 = nodeTagss[0][k] - 1;
		id2 = nodeTagss[0][k + 1] - 1;
		double xx = (Nodes[id1].getX() + Nodes[id2].getX()) / 2.0;
		double yy = (Nodes[id1].getY() + Nodes[id2].getY()) / 2.0;
		double zz = (Nodes[id1].getZ() + Nodes[id2].getZ()) / 2.0;
		gmsh::model::mesh::getElementsByCoordinates(xx, yy, zz, elementTag, dim2, strict);
		for (unsigned int kk = 0; kk < elementTag.size(); kk++)
		{
			//boundaryElementsIds.push_back(elementTag[kk] - 1);

			bottom_bound_egdes.push_back(boundary_edge(elementTag[kk] - 1, pair<int, int>{id1, id2}));
		}

	}
	allBoundEdges.push_back(bottom_bound_egdes);

	//right boundary
	dim = 1;
	dim2 = 2;
	strict = true;
	tag = l2;
	//vector<size_t> elementTag;
	vector<boundary_edge> right_bound_egdes;
	gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTagss, dim, tag);
	for (unsigned int k = 0; k < nodeTagss[0].size(); k += 2)
	{
		id1 = nodeTagss[0][k] - 1;
		id2 = nodeTagss[0][k + 1] - 1;
		double xx = (Nodes[id1].getX() + Nodes[id2].getX()) / 2.0;
		double yy = (Nodes[id1].getY() + Nodes[id2].getY()) / 2.0;
		double zz = (Nodes[id1].getZ() + Nodes[id2].getZ()) / 2.0;
		gmsh::model::mesh::getElementsByCoordinates(xx, yy, zz, elementTag, dim2, strict);
		for (unsigned int kk = 0; kk < elementTag.size(); kk++)
		{
			//boundaryElementsIds.push_back(elementTag[kk] - 1);

			right_bound_egdes.push_back(boundary_edge(elementTag[kk] - 1, pair<int, int>{id1, id2}));
		}

	}
	allBoundEdges.push_back(right_bound_egdes);


	//up boundary
	dim = 1;
	dim2 = 2;
	strict = true;
	tag = l3;
	//vector<size_t> elementTag;
	vector<boundary_edge> up_bound_egdes;
	gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTagss, dim, tag);
	for (unsigned int k = 0; k < nodeTagss[0].size(); k += 2)
	{
		id1 = nodeTagss[0][k] - 1;
		id2 = nodeTagss[0][k + 1] - 1;
		double xx = (Nodes[id1].getX() + Nodes[id2].getX()) / 2.0;
		double yy = (Nodes[id1].getY() + Nodes[id2].getY()) / 2.0;
		double zz = (Nodes[id1].getZ() + Nodes[id2].getZ()) / 2.0;
		gmsh::model::mesh::getElementsByCoordinates(xx, yy, zz, elementTag, dim2, strict);
		for (unsigned int kk = 0; kk < elementTag.size(); kk++)
		{
			//boundaryElementsIds.push_back(elementTag[kk] - 1);

			up_bound_egdes.push_back(boundary_edge(elementTag[kk] - 1, pair<int, int>{id1, id2}));
		}

	}
	allBoundEdges.push_back(up_bound_egdes);

	//left boundary
	dim2 = 2;
	strict = true;
	tag = l4;
	//vector<size_t> elementTag;
	vector<boundary_edge> left_bound_egdes;
	gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTagss, dim, tag);
	for (unsigned int k = 0; k < nodeTagss[0].size(); k += 2)
	{
		id1 = nodeTagss[0][k] - 1;
		id2 = nodeTagss[0][k + 1] - 1;
		double xx = (Nodes[id1].getX() + Nodes[id2].getX()) / 2.0;
		double yy = (Nodes[id1].getY() + Nodes[id2].getY()) / 2.0;
		double zz = (Nodes[id1].getZ() + Nodes[id2].getZ()) / 2.0;
		gmsh::model::mesh::getElementsByCoordinates(xx, yy, zz, elementTag, dim2, strict);
		for (unsigned int kk = 0; kk < elementTag.size(); kk++)
		{
			//boundaryElementsIds.push_back(elementTag[kk] - 1);
			left_bound_egdes.push_back(boundary_edge(elementTag[kk] - 1, pair<int, int>{id1, id2}));
		}

	}
	allBoundEdges.push_back(left_bound_egdes);

	gmsh::finalize();

	//export_to_tecplot(Nodes, triangles, "mesh.plt");
	return true;
}

//assign the boundary value on nodes
//size of the nodes is 4
//nodes[0]:left boundary (x=-1)
//nodes[1]:bottom boundary (y=-1)
//nodes[2]:right boundary (x=1)
//nodes[3]:up boundary (y=1)
void assign_boundarynodes_value(vector<vector<boundary_node>> &nodes, vector<Point>Pb_test)
{
	if (ex_index == 1)
	{
		//example 1 in chapter 3
		//left boundary
		for (unsigned int i = 0; i < nodes[0].size(); ++i)
		{
			nodes[0][i].type = Dirichlet;
			nodes[0][i].value = vector<double>{ left_bound_value(Pb_test[nodes[0][i].glo_node_index]) };
		}

		//bottom boundary
		for (unsigned int i = 0; i < nodes[1].size(); ++i)
		{
			nodes[1][i].type = Dirichlet;
			nodes[1][i].value = vector<double>{ bottom_bound_value(Pb_test[nodes[1][i].glo_node_index]) };
		}
		//right boundary
		for (unsigned int i = 0; i < nodes[2].size(); ++i)
		{
			nodes[2][i].type = Dirichlet;
			nodes[2][i].value = vector<double>{ right_bound_value(Pb_test[nodes[2][i].glo_node_index]) };
		}
		//up boundary
		for (unsigned int i = 0; i < nodes[3].size(); ++i)
		{
			nodes[3][i].type = Dirichlet;
			nodes[3][i].value = vector<double>{ up_bound_value(Pb_test[nodes[3][i].glo_node_index]) };
		}

	}

	else if (ex_index == 2 || 3 == ex_index)
	{
		//example 2 in chapter 3
		for (unsigned int i = 0; i < nodes[0].size(); ++i)
		{
			nodes[0][i].type = Dirichlet;
			nodes[0][i].value = vector<double>{ left_bound_value(Pb_test[nodes[0][i].glo_node_index]) };
		}

		for (unsigned int i = 0; i < nodes[1].size(); ++i)
		{
			nodes[1][i].type = Dirichlet;
			nodes[1][i].value = vector<double>{ right_bound_value(Pb_test[nodes[1][i].glo_node_index]) };
		}

		for (unsigned int i = 0; i < nodes[2].size(); ++i)
		{
			nodes[2][i].type = Dirichlet;
			nodes[2][i].value = vector<double>{ up_bound_value(Pb_test[nodes[2][i].glo_node_index]) };
		}
	}
	



}

/*
int dumpLine(FILE* fn)
{
	int c;
	int n = 0;//n=0 means the start line
	while ((c = fgetc(fn)) != '\n')
	{
		if (c == EOF)
		{
			n = 1;//the end of the file
			break;
		}
	}
	return n;
}*/

/*
bool read_mesh_from_file(string filename, vector<Point> &Nodes, vector<Triangle> &Triangles, vector<vector<boundary_node>> &allBoundNodes,
	vector<vector<boundary_edge>> &allBoundEdges, ELE_TYPE &ele_type)
{
	string suffixStr = filename.substr(filename.find_last_of('.') + 1);//get the suffix of the mesh file
	if (suffixStr == "neu")
	{
		int node_number;
		int ele_number;

		FILE *fn;
		if (fopen_s(&fn, filename.c_str(), "r"))
		{
			cout << "\n error: open mesh file error\n";
			return false;
		}
		else
		{
			try
			{
				dumpLine(fn);
				dumpLine(fn);
				dumpLine(fn);
				dumpLine(fn);
				dumpLine(fn);
				dumpLine(fn);
				int temp1, temp2, temp3, temp4;
				fscanf_s(fn, "%d%d%d%d%d%d\n", &node_number, &ele_number, &temp1, &temp2, &temp3, &temp4);
				dumpLine(fn);
				dumpLine(fn);

				Nodes.resize(node_number);
				Triangles.resize(ele_number);

				//	double t1;
				for (int i = 0; i < node_number; ++i)
				{
					int index;
					double x, y;
					fscanf_s(fn, "%d %lf %lf \n", &index, &x, &y);//read the node coordinates
					Nodes[i] = Point(x, y, 0);
				}
				dumpLine(fn);
				dumpLine(fn);

				for (int i = 0; i < ele_number; ++i)
				{
					char str_line[100];
					fgets(str_line, 100, fn);

					const char *delim = ": ";
					char *tempChar;
					vector<int> p;
					p.push_back(atoi(strtok(str_line, delim)));
					while (tempChar = strtok(NULL, delim))
					{
						p.push_back(atoi(tempChar));
					}
					ele_type = ELE_TYPE(p[2]);//单元类型

					for (int j = 3; j < p.size(); ++j)
					{
						Triangles[i].append(p[j] - 1);
					}

				}

				dumpLine(fn);
				//dumpLine(fn);
				char str1[30];
				int nTempBoundNum;
				char charTemp[10];
				char charTepm1[10];
				int GroupNum = 0;
				while (fscanf(fn, "%s %s %s\n", charTemp, str1, charTepm1))
				{
					string str = str1;
					if (str == "GROUP")
					{
						fscanf(fn, "%s%d%s%d%s%d%s%d\n", str1, &temp1, charTepm1, &nTempBoundNum, charTemp, &temp2, charTemp, &temp3);
						dumpLine(fn);
						dumpLine(fn);
						int p[10];
						for (int i = 0; i < nTempBoundNum / 10; ++i)
						{
							fscanf_s(fn, "%d%d%d%d%d%d%d%d%d%d", &p[0],
								&p[1],
								&p[2],
								&p[3],
								&p[4],
								&p[5],
								&p[6],
								&p[7],
								&p[8],
								&p[9]);

						}
						int nTempGroupRow = nTempBoundNum - nTempBoundNum % 10;
						for (int i = 0; i < nTempBoundNum % 10; ++i)
						{
							int temp;
							fscanf_s(fn, "%d", &temp);
						}
						dumpLine(fn);
						dumpLine(fn);
						GroupNum++;

					}
					else
					{
						break;
					}
				}

				int ntemp = 0;
				int data_type;//0=node,1=element
				while (fscanf(fn, "%s%d%d%d%d\n", str1, &data_type, &nTempBoundNum, &temp2, &temp3))
				{
					vector<boundary_node> boundarynodes;
					vector<boundary_edge> boundaryedges;
					if (data_type == 0)//the boundary is recorded by nodes
					{
						for (int i = 0; i < nTempBoundNum; ++i)
						{
							boundary_node bound_node;
							int index;
							fscanf_s(fn, "%d\n", &index);
							bound_node.glo_node_index = index - 1;
							boundarynodes.push_back(bound_node);
						}
					}
					else if (data_type == 1)//the boundary is recorded by element
					{
						for (int i = 0; i < nTempBoundNum; ++i)
						{
							boundary_edge boun_edge;
							int index;
							fscanf_s(fn, "%d%d%d\n", &index, &temp1, &temp2);//读取上边界面上的单元编号
							boun_edge.ele_index = index - 1;
							switch (temp2)
							{
							case 1:
								if (ELE_3 == ele_type)
								{
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[1]);
								}
								else if (ELE_6 == ele_type)
								{
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[1]);
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[2]);
								}

								break;
							case 2:
								if (ELE_3 == ele_type)
								{
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[1]);
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[2]);
								}
								else if (ELE_6 == ele_type)
								{
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[3]);
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[4]);
								}
								break;
							case 3:
								if (ELE_3 == ele_type)
								{
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[0]);
								}
								else if (ELE_6 == ele_type)
								{
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[4]);
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[5]);
									boun_edge.p.push_back(Triangles[boun_edge.ele_index].getIds()[0]);
								}
								break;
							default:
								break;
							}

						}
					}
					allBoundNodes.push_back(boundarynodes);
					allBoundEdges.push_back(boundaryedges);

					ntemp += 1;
					dumpLine(fn);
					if (dumpLine(fn))
						break;
				}

				fclose(fn);

				return true;
			}
			catch (exception& e)
			{
				cout << e.what() << endl;
			}


		}
	}
	else if (suffixStr == "msh")
	{

	}
}*/

//read mesh from file
//support the .neu mesh generated by Gambit and .msh file generated by gmsh
//return:
//mesh nodes and triangles,boundary nodes and boundary edges
bool read_mesh_from_file(string filename, vector<Point> &Nodes, vector<Triangle> &Triangles, vector<vector<boundary_node>> &allBoundNodes,
	vector<vector<boundary_edge>> &allBoundEdges, ELE_TYPE &ele_type)
{
	string suffixStr = filename.substr(filename.find_last_of('.') + 1);//get the suffix of the mesh file
	if (suffixStr == "neu")
	{
		int node_number;
		int ele_number;

		ifstream f(filename, ios::in);;
		if (!f.is_open())
		{
			cout << "\n error: open mesh file error\n";
			return false;
		}
		else
		{
			try
			{
				string s = "";
				for (int i = 0; i < 6; ++i)//skip 6 line for the file header
				{
					getline(f, s);
				}
				int NGRPS, NBSETS, NDFCD, NDFVL;
				f >> node_number >> ele_number >> NGRPS >> NBSETS >> NDFCD >> NDFVL;

				for (int i = 0; i < 3; ++i)//As >> operator does not handle the '\r' in the last line
				{
					getline(f, s);
				}

				Nodes.resize(node_number);
				Triangles.resize(ele_number);

				//	double t1;
				for (int i = 0; i < node_number; ++i)
				{
					int index;
					double x, y;
					f >> index >> x >> y;//read the node coordinates
					Nodes[i] = Point(x, y, 0);
				}
				for (int i = 0; i < 3; ++i)
				{
					getline(f, s);
				}

				for (int i = 0; i < ele_number; ++i)
				{
					int id; // no check is done to see that they are monotonously increasing
					int ntype, ndp;
					f >> id >> ntype >> ndp;
					ele_type = ELE_TYPE(ndp);
					Triangles[i].set_ele_type(ELE_TYPE(ndp));//set the element type, it is also the basis type
					for (int j = 0; j < ndp; ++j)
					{
						int list;
						f >> list;
						Triangles[i].append(list - 1);
					}
				}

				getline(f, s);
				getline(f, s);


				char str1[30];
				int nTempBoundNum;
				char charTemp[10];
				char charTepm1[10];
				int GroupNum = 0;
				int temp1, temp2, temp3;
				for (int i = 0; i < NGRPS; ++i)
				{
					f >> charTemp >> str1 >> charTepm1;
					string str = str1;
					if (str == "GROUP")
					{
						f >> str1 >> temp1 >> charTepm1 >> nTempBoundNum >> charTemp >> temp2 >> charTemp >> temp3;
						for (int i = 0; i < 3; ++i)
						{
							getline(f, s);
						}
						int p[10];
						for (int i = 0; i < nTempBoundNum / 10; ++i)
						{
							f >> p[0] >> p[1] >> p[2] >> p[3] >> p[4] >> p[5] >> p[6] >> p[7] >> p[8] >> p[9];

						}
						int nTempGroupRow = nTempBoundNum - nTempBoundNum % 10;
						for (int i = 0; i < nTempBoundNum % 10; ++i)
						{
							int temp;
							f >> temp;
						}
						for (int i = 0; i < 3; ++i)
						{
							getline(f, s);
						}
						GroupNum++;
					}
				}

				int ntemp = 0;
				int data_type;//0=node,1=element
				for (int i = 0; i < NBSETS; ++i)
				{
					f >> str1 >> data_type >> nTempBoundNum >> temp2 >> temp3;


					if (data_type == 0)//the boundary is recorded by nodes
					{
						vector<boundary_node> boundarynodes;
						for (int i = 0; i < nTempBoundNum; ++i)
						{
							boundary_node bound_node;
							int index;
							f >> index;
							bound_node.glo_node_index = index - 1;
							boundarynodes.push_back(bound_node);
						}
						allBoundNodes.push_back(boundarynodes);
					}
					else if (data_type == 1)//the boundary is recorded by element
					{
						vector<boundary_edge> boundaryedges;
						for (int i = 0; i < nTempBoundNum; ++i)
						{
							boundary_edge boun_edge;
							int index;
							f >> index >> temp1 >> temp2;//read element node on the boundary
							boun_edge.ele_index = index - 1;
							switch (temp2)
							{
							case 1:
								if (ELE_3 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[1]);
								}
								else if (ELE_6 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[2]);
								}
								else if (ELE_4 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[1]);
								}
								else if (ELE_9 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[0]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[2]);
								}

								break;
							case 2:
								if (ELE_3 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[1]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[2]);
								}
								else if (ELE_6 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[4]);
								}
								else if (ELE_4 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[1]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[2]);
								}
								else if (ELE_9 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[4]);
								}
								break;
							case 3:
								if (ELE_3 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[0]);
								}
								else if (ELE_6 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[4]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[0]);
								}
								else if (ELE_4 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[2]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[3]);
								}
								else if (ELE_9 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[4]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[6]);
								}
								break;
							case 4://case 4 only exists when the element has four edges
								
								if (ELE_4 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[3]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[0]);
								}
								else if (ELE_9 == ele_type)
								{
									boun_edge.p.first = (Triangles[boun_edge.ele_index].getIds()[6]);
									boun_edge.p.second = (Triangles[boun_edge.ele_index].getIds()[0]);
								}
								break;
							default:
								break;
							}
							boundaryedges.push_back(boun_edge);
						}
						allBoundEdges.push_back(boundaryedges);
					}
					ntemp += 1;
					for (int i = 0; i < 3; ++i)
					{
						getline(f, s);
					}
				}

				f.close();

				return true;
			}
			catch (exception& e)
			{
				cout << e.what() << endl;
				return false;
			}


		}
	}
	else if (suffixStr == "msh")
	{

	}
	else 
	{
		cout << "\nUnknown mesh file" << endl;
		return false;
	}
	return false;

}


//assign the boundary value on edges 	
void assign_boundaryedges_value(vector<vector<boundary_edge>> &edges, vector<Point>Pb_test)
{
	//example 2 in chapter 3
	if (ex_index == 2)
	{
		for (unsigned int i = 0; i < edges[0].size(); ++i)
		{
			edges[0][i].type = Neumann;
			edges[0][i].value = bottom_bound_value;//the value
		}
	}
	else if (ex_index==3)
	{
		for (unsigned int i = 0; i < edges[0].size(); ++i)
		{
			edges[0][i].type = Robin;
			edges[0][i].value = bottom_bound_value;//the q value
			edges[0][i].r = r;//the r value
		}
	}

}

double FE_reference_basis_func_2D(Point xh, BASIS_TYPE basis_type, int basis_index, int basis_der_x, int basis_der_y)
{
	double x = xh.getX();
	double y = xh.getY();
	if (basis_type == TWO_DIM_LINEAR)
	{
		if (basis_index == 0)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return -xh.getX() - xh.getY() + 1;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return -1;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return -1;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return 0;
			}
			else
				return 0;

		}
		if (basis_index == 1)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return xh.getX();
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return 0;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return 1;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return 0;
			}
			else
				return 0;

		}
		if (basis_index == 2)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return xh.getY();
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return 1;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return 0;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return 0;
			}
			else
				return 0;

		}
	}
	else if (basis_type == TWO_DIM_QUADRATIC)
	{
		//It should be noted that the node order of gambit is different from the index of the local basis function given by Prof. He
		//The basis function index is defined by Prof. He as follow:
		//     3 
		//    / \
		//   6   5  
		//  /     \
		// 1___4___2
		//However, the order of the node in a element generated by gambit is defined as follow:
		//     5 
		//    / \
		//   6   4  
		//  /     \
		// 1___2___3
		if (basis_index == 0)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return 2 * pow(xh.getX(), 2) + 2 * pow(xh.getY(), 2) + 4 * xh.getX()*xh.getY() - 3 * xh.getY() - 3 * xh.getX() + 1;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return 4 * xh.getX() + 4 * xh.getY() - 3;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return 4 * xh.getX() + 4 * xh.getY() - 3;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return 4;
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return 4;
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return 4;
			}
			else
				return 0;

		}
		if (basis_index == 2)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return 2 * pow(xh.getX(), 2) - xh.getX();
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return 0;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return 4 * xh.getX() - 1;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return 0;
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return 0;
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return 4;
			}
			else
				return 0;

		}
		if (basis_index == 4)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return 2 * pow(xh.getY(), 2) - xh.getY();
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return 4 * xh.getY() - 1;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return 0;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return 0;
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return 4;
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return 0;
			}
			else
				return 0;

		}
		if (basis_index == 1)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return -4 * pow(xh.getX(), 2) - 4 * xh.getX()*xh.getY() + 4 * xh.getX();
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return -4 * xh.getX();
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return -8 * xh.getX() - 4 * xh.getY() + 4;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return -4;
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return 0;
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return -8;
			}
			else
				return 0;

		}
		if (basis_index == 3)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return 4 * xh.getX()*xh.getY();
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return 4 * xh.getX();
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return 4 * xh.getY();
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return 4;
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return 0;
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return 0;
			}
			else
				return 0;

		}

		if (basis_index == 5)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return -4 * pow(xh.getY(), 2) - 4 * xh.getX()*xh.getY() + 4 * xh.getY();
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return -4 * xh.getX() - 8 * xh.getY() + 4;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return -4 * xh.getY();
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return -4;
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return -8;
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return 0;
			}
			else
				return 0;

		}
	}
	else if (basis_type == TWO_DIM_BILINEAR)
	{

		//define basis function on reference element([-1,-1],[1,-1],[1,1],[-1,1])

		if (basis_index == 0)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return (1 - xh.getX() - xh.getY() + xh.getX()*xh.getY()) / 4.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return (-1 + xh.getY()) / 4.0;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return (-1 + xh.getX()) / 4.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return 1.0 / 4.0;
			}
			else
				return 0;
		}
		if (basis_index == 1)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return (1 + xh.getX() - xh.getY() - xh.getX()*xh.getY()) / 4.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return (1 - xh.getY()) / 4.0;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return (-1 - xh.getX()) / 4.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return -1.0 / 4.0;
			}

			else
				return 0;
		}
		if (basis_index == 2)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return (1 + xh.getX() + xh.getY() + xh.getX()*xh.getY()) / 4.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return (1 + xh.getY()) / 4.0;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return (1 + xh.getX()) / 4.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return 1.0 / 4.0;
			}
			else
				return 0;
		}
		if (basis_index == 3)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return (1 - xh.getX() + xh.getY() - xh.getX()*xh.getY()) / 4.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return (-1 - xh.getY()) / 4.0;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return (1 - xh.getX()) / 4.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return -1.0 / 4.0;
			}
			else
				return 0;
		}
	}
	else if (basis_type == TWO_DIM_BIQUADRATIC)
	{
		//The reference basis functions and their derivatives are obtained by the symbolic calculation with python in the file "basis_fun_biquadratic_2D.py"
		if (basis_index == 0)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return pow(xh.getX(), 2) * pow(xh.getY(), 2) / 4 - pow(xh.getX(), 2) * xh.getY() / 4 - xh.getX() * pow(xh.getY(), 2) / 4 + xh.getX() * xh.getY() / 4;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return xh.getX() *  pow(xh.getY(), 2) / 2 - xh.getX() * xh.getY() / 2 - pow(xh.getY(), 2) / 4 + xh.getY() / 4;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return pow(xh.getX(), 2) * xh.getY() / 2 - pow(xh.getX(), 2) / 4 - xh.getX() * xh.getY() / 2 + xh.getX() / 4;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return xh.getX() * xh.getY() - xh.getX() / 2 - xh.getY() / 2 + 1.0 / 4.0;
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return xh.getY() * (xh.getY() - 1) / 2;
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return xh.getX() * (xh.getX() - 1) / 2;
			}
			else if (basis_der_x == 3 && basis_der_y == 0)
			{
				return 0;
			}
			else if (basis_der_x == 0 && basis_der_y == 3)
			{
				return 0;
			}
			else if (basis_der_x == 2 && basis_der_y == 1)
			{
				return xh.getY() - 1 / 2;
			}
			else if (basis_der_x == 1 && basis_der_y == 2)
			{
				return xh.getX() - 1 / 2;
			}
			else if (basis_der_x == 2 && basis_der_y == 2)
			{
				return 1;
			}
			else
				return 0;
		}
		else if (basis_index == 1)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return -pow(xh.getX(), 2) * pow(xh.getY(), 2) / 2 + pow(xh.getX(), 2) * xh.getY() / 2 + pow(xh.getY(), 2) / 2 - xh.getY() / 2;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return -xh.getX() * pow(xh.getY(), 2) + xh.getX() * xh.getY();
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return -pow(xh.getX(), 2) * xh.getY() + pow(xh.getX(), 2) / 2 + xh.getY() - 1.0 / 2.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return xh.getX() * (1 - 2 * xh.getY());
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return xh.getY() * (1 - xh.getY());
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return 1 - pow(xh.getX(), 2);
			}
			else if (basis_der_x == 3 && basis_der_y == 0)
			{
				return 0;
			}
			else if (basis_der_x == 0 && basis_der_y == 3)
			{
				return 0;
			}
			else if (basis_der_x == 2 && basis_der_y == 1)
			{
				return 1 - 2 * xh.getY();
			}
			else if (basis_der_x == 1 && basis_der_y == 2)
			{
				return -2 * xh.getY();
			}
			else if (basis_der_x == 2 && basis_der_y == 2)
			{
				return -2.0;
			}
			else
				return 0;
		}
		else if (basis_index == 2)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return pow(xh.getX(), 2) * pow(xh.getY(), 2) / 4 - pow(xh.getX(), 2) * xh.getY() / 4 + xh.getX() * pow(xh.getY(), 2) / 4 - xh.getX() * xh.getY() / 4;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return xh.getX() * pow(xh.getY(), 2) / 2 - xh.getX() * xh.getY() / 2 + pow(xh.getY(), 2) / 4 - xh.getY() / 4;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return pow(xh.getX(), 2) * xh.getY() / 2 - pow(xh.getX(), 2) / 4 + xh.getX() * xh.getY() / 2 - xh.getX() / 4;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return xh.getX() * xh.getY() - xh.getX() / 2 + xh.getY() / 2 - 1.0 / 4.0;
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return xh.getY() * (xh.getY() - 1) / 2;
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return xh.getX() * (xh.getX() + 1) / 2;
			}
			else if (basis_der_x == 3 && basis_der_y == 0)
			{
				return 0;
			}
			else if (basis_der_x == 0 && basis_der_y == 3)
			{
				return 0;
			}
			else if (basis_der_x == 2 && basis_der_y == 1)
			{
				return xh.getY() - 1.0 / 2.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 2)
			{
				return xh.getX() + 1.0 / 2.0;
			}
			else if (basis_der_x == 2 && basis_der_y == 2)
			{
				return 1;
			}
			else
				return 0;
		}
		else if (basis_index == 3)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return -pow(xh.getX(), 2) * pow(xh.getY(), 2) / 2 + pow(xh.getX(), 2) / 2 - xh.getX() * pow(xh.getY(), 2) / 2 + xh.getX() / 2;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return -xh.getX() * pow(xh.getY(), 2) + xh.getX() - pow(xh.getY(), 2) / 2 + 1.0 / 2.0;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return -pow(xh.getX(), 2)* xh.getY() - xh.getX() * xh.getY();
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return -xh.getY() * (2 * xh.getX() + 1);
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return 1 - pow(xh.getY(), 2);
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return -xh.getX() * (xh.getX() + 1);
			}
			else if (basis_der_x == 3 && basis_der_y == 0)
			{
				return 0;
			}
			else if (basis_der_x == 0 && basis_der_y == 3)
			{
				return 0;
			}
			else if (basis_der_x == 2 && basis_der_y == 1)
			{
				return -2 * xh.getY();
			}
			else if (basis_der_x == 1 && basis_der_y == 2)
			{
				return -(2 * xh.getX() + 1.0);
			}
			else if (basis_der_x == 2 && basis_der_y == 2)
			{
				return -2;
			}
			else
				return 0;
		}
		else if (basis_index == 4)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return pow(x, 2) * pow(y, 2) / 4 + pow(x, 2) * y / 4 + x * pow(y, 2) / 4 + x * y / 4;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return x * pow(y, 2) / 2 + x * y / 2 + pow(y, 2) / 4 + y / 4;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return pow(x, 2) * y / 2 + pow(x, 2) / 4 + x * y / 2 + x / 4;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return x * y + x / 2 + y / 2 + 1.0 / 4.0;
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return y * (y + 1) / 2;
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return x * (x + 1) / 2;
			}
			else if (basis_der_x == 3 && basis_der_y == 0)
			{
				return 0;
			}
			else if (basis_der_x == 0 && basis_der_y == 3)
			{
				return 0;
			}
			else if (basis_der_x == 2 && basis_der_y == 1)
			{
				return y + 1 / 2.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 2)
			{
				return x + 1 / 2.0;
			}
			else if (basis_der_x == 2 && basis_der_y == 2)
			{
				return 1;
			}
			else
				return 0;
		}

		else if (basis_index == 5)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return -pow(x, 2)* pow(y, 2) / 2 - pow(x, 2) * y / 2 + pow(y, 2) / 2 + y / 2;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return -x * pow(y, 2) - x * y;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return -pow(x, 2) * y - pow(x, 2) / 2 + y + 1.0 / 2.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return -x * (2 * y + 1);
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return -y * (y + 1);
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return 1 - pow(x, 2);
			}
			else if (basis_der_x == 3 && basis_der_y == 0)
			{
				return 0;
			}
			else if (basis_der_x == 0 && basis_der_y == 3)
			{
				return 0;
			}
			else if (basis_der_x == 2 && basis_der_y == 1)
			{
				return -(2 * y + 1);
			}
			else if (basis_der_x == 1 && basis_der_y == 2)
			{
				return -2 * x;
			}
			else if (basis_der_x == 2 && basis_der_y == 2)
			{
				return -2;
			}
			else
				return 0;
		}
		else if (basis_index == 6)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return pow(x, 2) * pow(y, 2) / 4 + pow(x, 2) * y / 4 - x * pow(y, 2) / 4 - x * y / 4;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return x * pow(y, 2) / 2 + x * y / 2 - pow(y, 2) / 4 - y / 4;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return pow(x, 2) * y / 2 + pow(x, 2) / 4 - x * y / 2 - x / 4;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return x * y + x / 2 - y / 2 - 1 / 4;
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return y * (y + 1) / 2;
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return x * (x - 1) / 2;
			}
			else if (basis_der_x == 3 && basis_der_y == 0)
			{
				return 0;
			}
			else if (basis_der_x == 0 && basis_der_y == 3)
			{
				return 0;
			}
			else if (basis_der_x == 2 && basis_der_y == 1)
			{
				return y + 1.0 / 2.0;
			}
			else if (basis_der_x == 1 && basis_der_y == 2)
			{
				return x - 1.0 / 2.0;
			}
			else if (basis_der_x == 2 && basis_der_y == 2)
			{
				return 1;
			}
			else
				return 0;
		}
		else if (basis_index == 7)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return -pow(x, 2) * pow(y, 2) / 2 + pow(x, 2) / 2 + x * pow(y, 2) / 2 - x / 2;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return -x * pow(y, 2) + x + pow(y, 2) / 2 - 1.0 / 2.0;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return -pow(x, 2) * y + x * y;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return y * (1 - 2 * x);
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return 1 - pow(y, 2);
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return x * (1 - x);
			}
			else if (basis_der_x == 3 && basis_der_y == 0)
			{
				return 0;
			}
			else if (basis_der_x == 0 && basis_der_y == 3)
			{
				return 0;
			}
			else if (basis_der_x == 2 && basis_der_y == 1)
			{
				return -2 * y;
			}
			else if (basis_der_x == 1 && basis_der_y == 2)
			{
				return 1 - 2 * x;
			}
			else if (basis_der_x == 2 && basis_der_y == 2)
			{
				return -2;
			}
			else
				return 0;
		}
		else if (basis_index == 8)
		{
			if (basis_der_x == 0 && basis_der_y == 0)
			{
				return pow(x, 2) * pow(y, 2) - pow(x, 2) - pow(y, 2) + 1;
			}
			else if (basis_der_x == 1 && basis_der_y == 0)
			{
				return 2 * x*pow(y, 2) - 2 * x;
			}
			else if (basis_der_x == 0 && basis_der_y == 1)
			{
				return 2 * pow(x, 2) * y - 2 * y;
			}
			else if (basis_der_x == 1 && basis_der_y == 1)
			{
				return 4 * x*y;
			}
			else if (basis_der_x == 2 && basis_der_y == 0)
			{
				return 2 * (pow(y, 2) - 1);
			}
			else if (basis_der_x == 0 && basis_der_y == 2)
			{
				return 2 * (pow(x, 2) - 1);
			}
			else if (basis_der_x == 3 && basis_der_y == 0)
			{
				return 0;
			}
			else if (basis_der_x == 0 && basis_der_y == 3)
			{
				return 0;
			}
			else if (basis_der_x == 2 && basis_der_y == 1)
			{
				return 4 * y;
			}
			else if (basis_der_x == 1 && basis_der_y == 2)
			{
				return 4 * x;
			}
			else if (basis_der_x == 2 && basis_der_y == 2)
			{
				return 4;
			}
			else
				return 0;
		}

	}
	return 0;
}

//local basis function for one dimension
//basis_type: define the basis type
//vertices: the coordinates of the vertices
//basis_index: index of the basis function, first one and second one for 1D problem
//basis_der_x: derivative order of the basis function for x
//The definition of basis local function for both test and trial

double FE_basis_local_func_2D(Point x, vector<Point> vertices, BASIS_TYPE basis_type, int basis_index, int basis_der_x, int basis_der_y)
{
	if (TWO_DIM_LINEAR == basis_type || TWO_DIM_QUADRATIC == basis_type)
	{
		double x1 = vertices[0].getX();
		double y1 = vertices[0].getY();
		double x2 = vertices[1].getX();
		double y2 = vertices[1].getY();
		double x3 = vertices[2].getX();
		double y3 = vertices[2].getY();

		double J = fabs((x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1));
		Point xh;
		xh.setX(((y3 - y1)*(x.getX() - x1) - (x3 - x1)*(x.getY() - y1)) / J);
		xh.setY((-(y2 - y1)*(x.getX() - x1) + (x2 - x1)*(x.getY() - y1)) / J);
		xh.setZ(0);

		if (0 == basis_der_x && 0 == basis_der_y)//zero derivative of the basis function
		{
			return FE_reference_basis_func_2D(xh, basis_type, basis_index, basis_der_x, basis_der_y);
		}
		else if (1 == basis_der_x && 0 == basis_der_y)//first derivation of the basis function
		{
			return (y3 - y1) / J * FE_reference_basis_func_2D(xh, basis_type, basis_index, 1, 0)
				+ (y1 - y2) / J * FE_reference_basis_func_2D(xh, basis_type, basis_index, 0, 1);
		}
		else if (0 == basis_der_x && 1 == basis_der_y)
		{
			return  (x1 - x3) / J * FE_reference_basis_func_2D(xh, basis_type, basis_index, 1, 0)
				+ (x2 - x1) / J * FE_reference_basis_func_2D(xh, basis_type, basis_index, 0, 1);
		}
		else if (2 == basis_der_x && 0 == basis_der_y)
		{
			return pow(y3 - y1, 2) / pow(J, 2)*FE_reference_basis_func_2D(xh, basis_type, basis_index, 2, 0) +
				2 * (y3 - y1)*(y1 - y2) / pow(J, 2)*FE_reference_basis_func_2D(xh, basis_type, basis_index, 1, 1) +
				pow(y1 - y2, 2) / pow(J, 2)*FE_reference_basis_func_2D(xh, basis_type, basis_index, 0, 2);
		}
		else if (0 == basis_der_x && 2 == basis_der_y)
		{
			return pow(x1 - x3, 2) / pow(J, 2)*FE_reference_basis_func_2D(xh, basis_type, basis_index, 2, 0) +
				2 * (x1 - x3)*(x2 - x1) / pow(J, 2)*FE_reference_basis_func_2D(xh, basis_type, basis_index, 1, 1) +
				pow(x2 - x1, 2) / pow(J, 2)*FE_reference_basis_func_2D(xh, basis_type, basis_index, 0, 2);
		}
		else if (1 == basis_der_x && 1 == basis_der_y)
		{
			return (x1 - x3) * (y3 - y1) / pow(J, 2)*FE_reference_basis_func_2D(xh, basis_type, basis_index, 2, 0) +
				(x1 - x3) * (y1 - y2) / pow(J, 2)*FE_reference_basis_func_2D(xh, basis_type, basis_index, 1, 1) +
				(x2 - x1) * (y3 - y1) / pow(J, 2)*FE_reference_basis_func_2D(xh, basis_type, basis_index, 1, 1) +
				(x2 - x1) * (y1 - y2) / pow(J, 2)*FE_reference_basis_func_2D(xh, basis_type, basis_index, 0, 2);
		}
		else if (basis_der_x > 2 || basis_der_y > 2)
		{
			return 0;
		}

		else
		{
			cout << "\n error: wrong input for basis derivative\n";
			return 0;
		}

	}
	else if (TWO_DIM_BILINEAR == basis_type || TWO_DIM_BIQUADRATIC == basis_type)
	{

		//The reference basis functions defined in the slides of Prof. He require that 
		//the A1,A2,A3,A4 are the left-lower,right-lower,right-upper, and left-upper vertices respectively
		//However, the mesh generated by the gambit may not meet the requirement.
		//So we developed a general affine transformation
		double x1 = vertices[0].getX();
		double y1 = vertices[0].getY();
		double x2 = vertices[1].getX();
		double y2 = vertices[1].getY();
		double x3 = vertices[2].getX();
		double y3 = vertices[2].getY();
		double x4 = vertices[3].getX();
		double y4 = vertices[3].getY();

		double a = (-2 * y2 + 2 * y3) / (x1*y2 - x1 * y3 - x2 * y1 + x2 * y3 + x3 * y1 - x3 * y2);
		double b = (2 * x2 - 2 * x3) / (x1*y2 - x1 * y3 - x2 * y1 + x2 * y3 + x3 * y1 - x3 * y2);
		double c = (x1*y2 - x1 * y3 - x2 * y1 - x2 * y3 + x3 * y1 + x3 * y2) /
			(x1*y2 - x1 * y3 - x2 * y1 + x2 * y3 + x3 * y1 - x3 * y2);
		double d = (2 * y1 - 2 * y2) / (x1*y2 - x1 * y3 - x2 * y1 + x2 * y3 + x3 * y1 - x3 * y2);
		double e = (-2 * x1 + 2 * x2) / (x1*y2 - x1 * y3 - x2 * y1 + x2 * y3 + x3 * y1 - x3 * y2);
		double f = (x1*y2 + x1 * y3 - x2 * y1 - x2 * y3 - x3 * y1 + x3 * y2) /
			(x1*y2 - x1 * y3 - x2 * y1 + x2 * y3 + x3 * y1 - x3 * y2);
		double J = fabs((a*e - b * d) / pow(a*e - b * d, 2));

		Point xh;
		xh.setX(a * x.getX() + b * x.getY() + c);
		xh.setY(d * x.getX() + e * x.getY() + f);
		xh.setZ(0);


		if (0 == basis_der_x && 0 == basis_der_y)//zero derivative of the basis function
		{
			return FE_reference_basis_func_2D(xh, basis_type, basis_index, basis_der_x, basis_der_y);
		}
		else if (1 == basis_der_x && 0 == basis_der_y)//first derivation of the basis function
		{
			return a * FE_reference_basis_func_2D(xh, basis_type, basis_index, 1, 0) +
				d * FE_reference_basis_func_2D(xh, basis_type, basis_index, 0, 1);
		}
		else if (0 == basis_der_x && 1 == basis_der_y)
		{
			return  b * FE_reference_basis_func_2D(xh, basis_type, basis_index, 1, 0) +
				e * FE_reference_basis_func_2D(xh, basis_type, basis_index, 0, 1);

		}
		else if (1 == basis_der_x && 1 == basis_der_y)
		{
			return
				a * b*FE_reference_basis_func_2D(xh, basis_type, basis_index, 2, 0) +
				b * d*FE_reference_basis_func_2D(xh, basis_type, basis_index, 1, 1) +
				a * e*FE_reference_basis_func_2D(xh, basis_type, basis_index, 1, 1) +
				d * e*FE_reference_basis_func_2D(xh, basis_type, basis_index, 0, 2);

		}


	}

	return 0;
}

double Gauss_quadrature_2D_trial_test(double(*coe)(Point), vector<double> Gauss_weights, vector<Point> Gauss_nodes, vector<Point> vertices,
	BASIS_TYPE basis_type_trial, int basis_index_trial, int basis_der_x_trial, int basis_der_y_trial,
	BASIS_TYPE basis_type_test, int basis_index_test, int basis_der_x_test, int basis_der_y_test)
{
	int Gpn = Gauss_weights.size();
	double int_value = 0;

	for (int k = 0; k < Gpn; ++k)
	{
		int_value += Gauss_weights[k] * coe(Gauss_nodes[k])*FE_basis_local_func_2D(Gauss_nodes[k], vertices, basis_type_trial, basis_index_trial, basis_der_x_trial, basis_der_y_trial)*FE_basis_local_func_2D(Gauss_nodes[k], vertices, basis_type_test, basis_index_test, basis_der_x_test, basis_der_y_test);
	}
	return int_value;
}
double Gauss_quadrature_2D_load_trial(double(*load)(Point), vector<double> Gauss_weights, vector<Point> Gauss_nodes, vector<Point> vertices,
	BASIS_TYPE basis_type_test, int basis_index_test, int basis_der_x_test, int basis_der_y_test)
{
	int Gpn = Gauss_weights.size();
	double int_value = 0;
	for (int k = 0; k < Gpn; ++k)
	{
		int_value += Gauss_weights[k] * load(Gauss_nodes[k])*FE_basis_local_func_2D(Gauss_nodes[k], vertices, basis_type_test, basis_index_test, basis_der_x_test, basis_der_y_test);
	}
	return int_value;
}

//generate Gauss points on a rectangle([-1,-1],[1,-1],[1,1],[-1,1])
void generate_Gauss_local_2D_Rec(vector<Point> vertices, int Gauss_type, vector<Point> &Gauss_nodes, vector<double> &Gauss_weights)
{
	//Gauss quadrature points in the reference rectangle
	//Refer to https://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules/quadrature_rules.html to get other type of numerical integration formula
	vector<double> w;
	vector<Point> p;
	if (9 == Gauss_type)
	{
		/*w = vector<double>{
			0.1111111111111111,
			0.4444444444444443,
			0.1111111111111111,
			0.4444444444444443,
			1.777777777777777,
			0.4444444444444443,
			0.1111111111111111,
			0.4444444444444443,
			0.1111111111111111 };
		p = vector<Point>{
			Point(-1.0,-1.0,0.0),
			Point(-1.0,0.0,0.0),
			Point(-1.0,1.0,0.0),
			Point(0.0,-1.0,0.0),
			Point(0.0,0.0,0.0),
			Point(0.0,1.0,0.0),
			Point(1.0,-1.0,0.0),
			Point(1.0,0.0,0.0),
			Point(1.0,1.0,0.0)};*/

			/*w = {
				.3086419753086420,
				.4938271604938272,
				.3086419753086420,
				.4938271604938272,
				.7901234567901235,
				.4938271604938272,
				.3086419753086420,
				.4938271604938272,
				.3086419753086420 };
			p = vector<Point>{
				Point(-.7745966692414835 -.7745966692414835,0),
				Point(-.7745966692414835,  .0000000000000000,0),
				Point(-.7745966692414835,  .7745966692414835,0),
				Point(.0000000000000000, -.7745966692414835,0),
				Point(.0000000000000000,  .0000000000000000,0),
				Point(.0000000000000000, .7745966692414835,0),
				Point(.7745966692414835, -.7745966692414835,0),
				Point(.7745966692414835, .0000000000000000,0),
				Point(.7745966692414835, .7745966692414835,0) };*/

		double gpt = pow(3.0 / 5.0, 1.0 / 2.0);
		double A1 = 25.0 / 81.0;
		double A2 = 40.0 / 81;
		double A3 = 64.0 / 81;
		w = {
			A1,
			A1,
			A1,
			A1,
			A2,
			A2,
			A2,
			A2,
			A3
		};
		p = vector<Point>{
			Point(-gpt,-gpt,0),
			Point(gpt, -gpt,0),
			Point(gpt,  gpt,0),
			Point(-gpt, gpt,0),
			Point(0,  -gpt,0),
			Point(gpt, 0,0),
			Point(0, gpt,0),
			Point(-gpt, 0,0),
			Point(0, 0,0) };

	}
	else
	{
		cout << "\nerror: Gauss quadrature type error\n";
	}

	//generate the Gauss points in the actual rectangle
	//The reference basis functions defined in the slides of Prof. He require that 
	//the A1,A2,A3,A4 are the left-lower,right-lower,right-upper, and left-upper vertices respectively
	//However, the mesh generated by the gambit may not meet the requirement.
	//So we developed a general affine transformation to transform any arbitrary rectangle to [-1,-1][1,1]
	double x1 = vertices[0].getX();
	double y1 = vertices[0].getY();
	double x2 = vertices[1].getX();
	double y2 = vertices[1].getY();
	double x3 = vertices[2].getX();
	double y3 = vertices[2].getY();
	double x4 = vertices[3].getX();
	double y4 = vertices[3].getY();

	double a = (-2 * y2 + 2 * y3) / (x1*y2 - x1 * y3 - x2 * y1 + x2 * y3 + x3 * y1 - x3 * y2);
	double b = (2 * x2 - 2 * x3) / (x1*y2 - x1 * y3 - x2 * y1 + x2 * y3 + x3 * y1 - x3 * y2);
	double c = (x1*y2 - x1 * y3 - x2 * y1 - x2 * y3 + x3 * y1 + x3 * y2) /
		(x1*y2 - x1 * y3 - x2 * y1 + x2 * y3 + x3 * y1 - x3 * y2);
	double d = (2 * y1 - 2 * y2) / (x1*y2 - x1 * y3 - x2 * y1 + x2 * y3 + x3 * y1 - x3 * y2);
	double e = (-2 * x1 + 2 * x2) / (x1*y2 - x1 * y3 - x2 * y1 + x2 * y3 + x3 * y1 - x3 * y2);
	double f = (x1*y2 + x1 * y3 - x2 * y1 - x2 * y3 - x3 * y1 + x3 * y2) /
		(x1*y2 - x1 * y3 - x2 * y1 + x2 * y3 + x3 * y1 - x3 * y2);
	double J = fabs(1 / (a*e - b * d));

	for (unsigned i = 0; i < w.size(); ++i)
	{
		Gauss_weights.push_back(w[i] * J);
		Gauss_nodes.push_back(
			Point(
			(b * f - b * p[i].getY() - c * e + e * p[i].getX()) / (a * e - b * d),
				(-a * f + a * p[i].getY() + c * d - d * p[i].getX()) / (a * e - b * d),
				0));
	}
}

//generate Gauss points on a triangle([0,0],[1,0],[0,1])
void generate_Gauss_local_2D_Tri(vector<Point> vertices, int Gauss_type, vector<Point> &Gauss_nodes, vector<double> &Gauss_weights)
{
	//Gauss quadrature points in the reference triangle
	vector<double> w;
	vector<Point> p;

	if (4 == Gauss_type)
	{
		w = vector<double>{ (1.0 - 1.0 / sqrt(3.0)) / 8.0,(1 - 1 / sqrt(3.0)) / 8,(1.0 + 1.0 / sqrt(3.0)) / 8,(1 + 1 / sqrt(3.0)) / 8.0 };
		p = vector<Point>{
			Point((1 / sqrt(3.0) + 1) / 2.0,(1.0 - 1.0 / sqrt(3.0))*(1.0 + 1 / sqrt(3.0)) / 4.0,0),
			Point((1 / sqrt(3.0) + 1) / 2.0,(1.0 - 1.0 / sqrt(3.0))*(1.0 - 1 / sqrt(3.0)) / 4.0,0),
			Point((-1 / sqrt(3.0) + 1) / 2.0,(1.0 + 1.0 / sqrt(3.0))*(1.0 + 1 / sqrt(3.0)) / 4.0,0),
			Point((-1 / sqrt(3.0) + 1) / 2.0,(1.0 + 1.0 / sqrt(3.0))*(1.0 - 1 / sqrt(3.0)) / 4.0,0)
		};
	}

	else if (9 == Gauss_type)
	{
		w = vector<double>{
			64.0 / 81.0 * (1.0 - 0.0) / 8,
			100.0 / 324.0 * (1.0 - sqrt(3.0 / 5.0)) / 8.0,
			100.0 / 324.0 * (1 - sqrt(3.0 / 5.0)) / 8.0,
			100.0 / 324.0 * (1 + sqrt(3.0 / 5.0)) / 8.0,
			100.0 / 324.0 * (1.0 + sqrt(3.0 / 5.0)) / 8,
			40.0 / 81.0 * (1.0 - 0.0) / 8.0,
			40.0 / 81.0 * (1.0 - 0.0) / 8.0,
			40.0 / 81.0 * (1.0 - sqrt(3.0 / 5.0)) / 8.0,
			40.0 / 81.0 * (1.0 + sqrt(3.0 / 5.0)) / 8.0 };
		p = vector<Point>{
			Point((1. + 0) / 2.,(1. - 0.)*(1. + 0.) / 4.,0),
			Point((1. + sqrt(3. / 5.)) / 2.,(1. - sqrt(3. / 5.))*(1. + sqrt(3. / 5.)) / 4.,0),
			Point((1. + sqrt(3. / 5.)) / 2,(1 - sqrt(3. / 5.))*(1. - sqrt(3. / 5.)) / 4.,0),
			Point((1. - sqrt(3. / 5.)) / 2,(1 + sqrt(3. / 5.))*(1 + sqrt(3. / 5.)) / 4.,0),
			Point((1. - sqrt(3. / 5.)) / 2,(1 + sqrt(3. / 5.))*(1. - sqrt(3. / 5.)) / 4.,0),
			Point((1. + 0) / 2.,(1. - 0)*(1. + sqrt(3. / 5.)) / 4.,0),
			Point((1. + 0) / 2.,(1. - 0)*(1. - sqrt(3. / 5.)) / 4.,0),
			Point((1. + sqrt(3. / 5.)) / 2,(1 - sqrt(3. / 5.))*(1. + 0) / 4.,0),
			Point((1. - sqrt(3. / 5.)) / 2,(1 + sqrt(3. / 5.))*(1. + 0) / 4.,0) };
	}
	else
	{
		cout << "\n error: Gauss quadrature type error\n";
		return;
	}

	//generate the Gauss points in the actual triangle
	//affine transformation
	double x1 = vertices[0].getX();
	double y1 = vertices[0].getY();
	double x2 = vertices[1].getX();
	double y2 = vertices[1].getY();
	double x3 = vertices[2].getX();
	double y3 = vertices[2].getY();

	double J = fabs((x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1));

	for (unsigned int i = 0; i < w.size(); ++i)
	{
		Gauss_weights.push_back(w[i] * J);
		Gauss_nodes.push_back(
			Point(
				x1 + (x2 - x1)*p[i].getX() + (x3 - x1)*p[i].getY(),
				y1 + (y2 - y1)*p[i].getX() + (y3 - y1)*p[i].getY(),
				0));
	}

}

SpMat assemble_matrix_2D(pair<int, int>matrix_size, vector<Point>P, vector<Triangle>T, vector<Triangle>Tb_trial, vector<Triangle>Tb_test,
	BASIS_TYPE basis_type_trial, int basis_der_x_trial, int basis_der_y_trial, BASIS_TYPE basis_type_test, int basis_der_x_test, int basis_der_y_test, int Gauss_type)
{
	int estimation_of_entries = Tb_trial[0].get_ele_node_num()*Tb_trial.size(); //estimation of the non-zero values
	vector<Tri> tripletList;//a triple list for the matrix to improve the efficiency of sparse matrix manipulation
	tripletList.reserve(estimation_of_entries);

	for (unsigned int n = 0; n < Tb_trial.size(); ++n)
	{
		vector<Point> vertices;//vertices (coordinates) of nth element(mesh nodes not finite element nodes)		
		for (int i = 0; i < T[n].get_ele_node_num(); ++i)
		{
			vertices.push_back(P[T[n].getIds()[i]]);
		}

		vector<double> Gauss_weights;
		vector<Point> Gauss_nodes;
		//generate the Gauss integration points
		//Gauss integration points over a triangle is different from a rectangle
		if (TWO_DIM_LINEAR == basis_type_test || TWO_DIM_QUADRATIC == basis_type_test)
		{
			generate_Gauss_local_2D_Tri(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		}
		else if (TWO_DIM_BILINEAR == basis_type_test || TWO_DIM_BIQUADRATIC == basis_type_test)
		{
			generate_Gauss_local_2D_Rec(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		}
		/*if (ELE_3 == Tb_trial[n].getIds().size() || ELE_6 == Tb_trial[n].getIds().size())
		{
			generate_Gauss_local_2D_Tri(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		}
		else if (ELE_4 == Tb_trial[n].getIds().size() || ELE_9 == Tb_trial[n].getIds().size())
		{
			generate_Gauss_local_2D_Rec(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		}*/

		for (unsigned int alpha = 0; alpha < Tb_trial[n].getIds().size(); ++alpha)
			for (unsigned int beta = 0; beta < Tb_test[n].getIds().size(); ++beta)
			{
				double int_value_x = Gauss_quadrature_2D_trial_test(coe_function, Gauss_weights, Gauss_nodes, vertices,
					basis_type_trial, alpha, basis_der_x_trial, 0,
					basis_type_test, beta, basis_der_x_test, 0);
				double ini_value_y = Gauss_quadrature_2D_trial_test(coe_function, Gauss_weights, Gauss_nodes, vertices,
					basis_type_trial, alpha, 0, basis_der_y_trial,
					basis_type_test, beta, 0, basis_der_y_test);

				int i = Tb_test[n].getIds()[beta];
				int j = Tb_trial[n].getIds()[alpha];
				tripletList.push_back(Tri(i, j, int_value_x + ini_value_y));
			}
	}
	SpMat A(matrix_size.first, matrix_size.second);
	A.setFromTriplets(tripletList.begin(), tripletList.end());
	return A;
}

Rsv assemble_vector_2D(pair<int, int>matrix_size, vector<Point>P, vector<Triangle>T, vector<Triangle>Tb_test, BASIS_TYPE basis_type_test, int basis_der_x_test, int basis_der_y_test, int Gauss_type)
{
	auto b = Rsv(matrix_size.first);

	for (int i = 0; i < matrix_size.first; ++i)
	{
		b(i) = 0;
	}

	for (unsigned int n = 0; n < Tb_test.size(); ++n)
	{
		vector<Point> vertices;//vertices of nth element
		vector<double> Gauss_weights;
		vector<Point> Gauss_nodes;

		for (int i = 0; i < T[n].get_ele_node_num(); ++i)
		{
			vertices.push_back(P[T[n].getIds()[i]]);
		}
		//generate the Gauss integration points
		if (TWO_DIM_LINEAR == basis_type_test || TWO_DIM_QUADRATIC == basis_type_test)
		{
			generate_Gauss_local_2D_Tri(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		}
		else if (TWO_DIM_BILINEAR == basis_type_test || TWO_DIM_BIQUADRATIC == basis_type_test)
		{
			generate_Gauss_local_2D_Rec(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		}

		/*if (ELE_3 == Tb_test[n].getIds().size() || ELE_6 == Tb_test[n].getIds().size())
		{
			generate_Gauss_local_2D_Tri(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		}
		else if (ELE_4 == Tb_test[n].getIds().size() || ELE_9 == Tb_test[n].getIds().size())
		{
			generate_Gauss_local_2D_Rec(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		}*/


		for (unsigned int beta = 0; beta < Tb_test[n].getIds().size(); ++beta)
		{
			double int_value = Gauss_quadrature_2D_load_trial(load_function, Gauss_weights, Gauss_nodes, vertices, basis_type_test, beta, basis_der_x_test, basis_der_y_test);

			b[Tb_test[n].getIds()[beta]] += int_value;
		}
	}
	return b;
}
void generate_Gauss_local_1D(vector<Point> vertices, int Gauss_type, vector<Point> &Gauss_nodes, vector<double> &Gauss_weights)
{
	//Gauss quadrature points on [-1,1]
	vector<double> w;
	vector<double> p;
	if (Gauss_type == 2)
	{
		w = vector<double>{ 1,1 };
		p = vector<double>{ -1 / sqrt(3.0), 1 / sqrt(3.0) };
	}

	else if (Gauss_type == 4)
	{
		w = vector<double>{ 0.3478548451, 0.3478548451, 0.6521451549, 0.6521451549 };
		p = vector<double>{ 0.8611363116, -0.8611363116, 0.3399810436, -0.3399810436 };
	}

	else if (Gauss_type == 8)
	{
		w = vector<double>{ 0.1012285363, 0.1012285363, 0.2223810345, 0.2223810345, 0.3137066459, 0.3137066459, 0.3626837834, 0.3626837834 };
		p = vector<double>{ 0.9602898565, -0.9602898565, 0.7966664774, -0.7966664774, 0.5255324099, -0.5255324099, 0.1834346425, -0.1834346425 };
	}
	else if (9==Gauss_type)
	{
		w = vector<double>{ 
			0.08127439, 0.18064816, 0.2606107,  0.31234708, 0.33023936,
			0.31234708,	0.2606107,  0.18064816, 0.08127439 };
		p = vector<double>{
			-0.96816024, -0.83603111, -0.61337143, -0.32425342,  0. ,
			0.32425342,	0.61337143,  0.83603111,  0.96816024 };
	}
	else
	{
		cout << "\n error: Gauss quadrature type error\n";
		return;
	}

	//Gauss quadrature point on [vertices[0],vertices[1]]
	for (unsigned int i = 0; i < w.size(); ++i)
	{

		Gauss_weights.push_back(sqrt(pow(vertices[0].getX() - vertices[1].getX(), 2) + pow(vertices[0].getY() - vertices[1].getY(), 2)) / 2.0 *w[i]);
		Gauss_nodes.push_back(
			Point(
			    (vertices[1].getX() - vertices[0].getX())*p[i] / 2 + (vertices[1].getX() + vertices[0].getX()) / 2.0,
				(vertices[1].getY() - vertices[0].getY())*p[i] / 2 + (vertices[1].getY() + vertices[0].getY()) / 2.0,
				0
			));
	}

}

double Gauss_quadrature_basis_1D_edge(double(*p_bar)(Point,double(*p)(Point)), double(*p)(Point),vector<double> Gauss_weights, vector<Point> Gauss_nodes, vector<Point> vertices, 
	BASIS_TYPE basis_type, int basis_index, int basis_der_x, int basis_der_y)
{
	int Gpn = Gauss_weights.size();
	double int_value = 0;
	for (int k = 0; k < Gpn; ++k)
	{
		int_value += Gauss_weights[k] * p_bar(Gauss_nodes[k],p)*FE_basis_local_func_2D(Gauss_nodes[k], vertices, basis_type, basis_index, basis_der_x, basis_der_y);
	}
	return int_value;
}

double Gauss_quadrature_basis_basis_1D_edge(double(*cr)(Point, double(*r)(Point)), double(*r)(Point), vector<double> Gauss_weights, vector<Point> Gauss_nodes, vector<Point> vertices, 
	BASIS_TYPE basis_type_trial, int basis_index_trial, int basis_der_x_trial, int basis_der_y_trial,
	BASIS_TYPE basis_type_test, int basis_index_test, int basis_der_x_test, int basis_der_y_test)
{
	int Gpn = Gauss_weights.size();
	double int_value = 0;
	for (int k = 0; k < Gpn; ++k)
	{
		int_value += Gauss_weights[k] * cr(Gauss_nodes[k], r)*FE_basis_local_func_2D(Gauss_nodes[k], vertices, basis_type_trial, basis_index_trial, basis_der_x_trial, basis_der_y_trial)
			*FE_basis_local_func_2D(Gauss_nodes[k], vertices, basis_type_test, basis_index_test, basis_der_x_test, basis_der_y_test);
	}
	return int_value;
}
double p_bar(Point x, double(*p)(Point))
{
	return coe_function(x)*p(x);
}

void treat_boundary(SpMat &A, Rsv &b, vector<vector<boundary_node>> boundarynodes, vector<vector<boundary_edge>> boundaryedges, vector<Point>P, vector<Triangle>T, vector<Point> Pb_test, vector<Triangle> Tb_test, BASIS_TYPE basis_type_test, vector<Triangle> Tb_trial, BASIS_TYPE basis_type_trial,int Gauss_type)
{
	for (unsigned int k = 0; k < boundarynodes.size(); ++k)
	{
		for (unsigned int l = 0; l < boundarynodes[k].size(); ++l)
		{
			if (Dirichlet == boundarynodes[k][l].type)
			{
				//set the A(i,i)=1 and A(i,j)=0, b=g(Pb(i)) to handle Dirichlet boundary
				//It is needed a loop over the column, and the manipulation of matrix is severely time consuming 
				//So we use the the second way to handle Dirichlet boundary

				/*int i = boundarynodes[k][l].glo_node_index;
				for (int j = 0; j < A.cols(); ++j)
				{
					A.insert(i, j) = 0;
				}
				A.insert(i, i) = 1;
				b[i] = boundarynodes[k][l].value[0];*/

				//another way to handle the Dirichlet boundary
				//b(i)=b=g(Pb(i))*A(i,i)*large_number
				//A(i,i)=A(i,i)*large_number
				//It is only needed to manipulate the diagonal element of the matrix

				int i = boundarynodes[k][l].glo_node_index;
				double large_number = 1e50;
				Eigen::SparseMatrix<double>::InnerIterator it(A, i);
				b(i) = it.value()*large_number*boundarynodes[k][l].value[0];
				A.coeffRef(i, i) = it.value()* large_number;
			}
		}

	}

	for (unsigned int k = 0; k < boundaryedges.size(); ++k)
	{
		for (unsigned int l = 0; l < boundaryedges[k].size(); ++l)
		{
			if (Neumann == boundaryedges[k][l].type)
			{
				int nk = boundaryedges[k][l].ele_index;
				int p1 = boundaryedges[k][l].p.first;
				int p2 = boundaryedges[k][l].p.second;

				//This vertices is the boundary edge's two nodes
				//It is for the generation of Gauss points on the edges
				vector<Point> vertices_ede{
					Pb_test[p1],
					Pb_test[p2]
				};
				vector<Point>Gauss_nodes;
				vector<double>Gauss_wights;
				generate_Gauss_local_1D(vertices_ede, Gauss_type, Gauss_nodes, Gauss_wights);

				//The basis function needs the vertices of the element not the edge
				//So it is needed to get the nk's nodes
				vector<Point> vertices;//vertices (coordinates) of nth element(mesh nodes not finite element nodes)		
				for (int i = 0; i < T[nk].get_ele_node_num(); ++i)
				{
					vertices.push_back(P[T[nk].getIds()[i]]);
				}
				for (unsigned int beta = 0; beta < Tb_test[nk].getIds().size(); ++beta)
				{
					double ini_value = Gauss_quadrature_basis_1D_edge(p_bar, boundaryedges[k][l].value,Gauss_wights, Gauss_nodes, vertices, basis_type_test, beta, 0, 0);
					b(Tb_test[nk].getIds()[beta]) += ini_value;
				}
			}
			else if (Robin == boundaryedges[k][l].type)
			{
				int nk = boundaryedges[k][l].ele_index;
				int p1 = boundaryedges[k][l].p.first;
				int p2 = boundaryedges[k][l].p.second;

				//This vertices is the boundary edge's two nodes
				//It is for the generation of Gauss points on the edges
				vector<Point> vertices_ede{
					Pb_test[p1],
					Pb_test[p2]
				};
				vector<Point>Gauss_nodes;
				vector<double>Gauss_wights;
				generate_Gauss_local_1D(vertices_ede, Gauss_type, Gauss_nodes, Gauss_wights);

				//The basis function needs the vertices of the element not the edge
				//So it is needed to get the nk's nodes
				vector<Point> vertices;//vertices (coordinates) of nth element(mesh nodes not finite element nodes)		
				for (int i = 0; i < T[nk].get_ele_node_num(); ++i)
				{
					vertices.push_back(P[T[nk].getIds()[i]]);
				}
				for (unsigned int beta = 0; beta < Tb_test[nk].getIds().size(); ++beta)
				{
					double ini_value = Gauss_quadrature_basis_1D_edge(p_bar, boundaryedges[k][l].value, Gauss_wights, Gauss_nodes, vertices, basis_type_test, beta, 0, 0);
					b(Tb_test[nk].getIds()[beta]) += ini_value;
					
					for (unsigned int alpha=0;alpha<Tb_trial[nk].getIds().size();++alpha)
					{
						double ini_value = Gauss_quadrature_basis_basis_1D_edge(p_bar, boundaryedges[k][l].r, Gauss_wights, Gauss_nodes, vertices,
							basis_type_trial, alpha, 0, 0, basis_type_test, beta, 0, 0);
						int i = Tb_test[nk].getIds()[beta];
						int j = Tb_trial[nk].getIds()[alpha];
						A.coeffRef(i, j) += ini_value;

					}
				}
			}
		}
		
	}
}
void set_basis_func_type(ELE_TYPE ele_type, BASIS_TYPE &basis_type_trial, BASIS_TYPE &basis_type_test)
{
	if (ELE_3 == ele_type)
	{
		basis_type_trial = TWO_DIM_LINEAR;
		basis_type_test = TWO_DIM_LINEAR;
	}
	else if (ELE_4 == ele_type)
	{

		basis_type_trial = TWO_DIM_BILINEAR;
		basis_type_test = TWO_DIM_BILINEAR;
	}
	else if (ELE_6 == ele_type)
	{
		basis_type_trial = TWO_DIM_QUADRATIC;
		basis_type_test = TWO_DIM_QUADRATIC;
	}
	else if (ELE_9 == ele_type)
	{
		basis_type_trial = TWO_DIM_BIQUADRATIC;
		basis_type_test = TWO_DIM_BIQUADRATIC;
	}
}
//integrate of the local solution
double local_FE_function_2D(vector<double>uh_vector, int number_of_local_func, Point x, vector<Point> vertices, BASIS_TYPE basis_type, int basis_der_x, int basis_der_y)
{
	double result = 0;
	for (int k = 0; k < number_of_local_func; ++k)
	{
		result += uh_vector[k] * FE_basis_local_func_2D(x, vertices, basis_type, k, basis_der_x, basis_der_y);

	}
	return result;
}

//Compute the error with Gauss integration
double Gauss_int_error_2D(double(*ana_sol)(Point, int, int), vector<double> uh_vector, vector<double> Gauss_weights, vector<Point> Gauss_nodes, vector<Point> vertices,
	BASIS_TYPE basis_type_trial, int der_x, int der_y, int number_of_local_func)
{
	int Gpn = Gauss_weights.size();
	double int_value_x = 0;
	double int_value_y = 0;

	for (int k = 0; k < Gpn; ++k)
	{
		int_value_x += Gauss_weights[k] * pow(ana_sol(Gauss_nodes[k], der_x, 0) - local_FE_function_2D(uh_vector, number_of_local_func, Gauss_nodes[k], vertices, basis_type_trial, der_x, 0), 2);
		int_value_y += Gauss_weights[k] * pow(ana_sol(Gauss_nodes[k], 0, der_y) - local_FE_function_2D(uh_vector, number_of_local_func, Gauss_nodes[k], vertices, basis_type_trial, 0, der_y), 2);

	}
	return int_value_x + int_value_y;
}


//Compute the L2 and H1 error
//der_x=0,der_y=0: L2
//der_x=1,der_y=1: H1
double compute_Hs_error(vector<Point> P, vector<Triangle>T, vector<Triangle>Tb_test,
	BASIS_TYPE basis_type_trial, int der_x, int der_y, int Gauss_type, vector<double> solution_vector)
{
	double error = 0;
	for (unsigned int n = 0; n < Tb_test.size(); ++n)
	{
		vector<Point> vertices;//vertices (coordinates) of nth element(mesh nodes not finite element nodes)		
		for (int i = 0; i < T[n].get_ele_node_num(); ++i)
		{
			vertices.push_back(P[T[n].getIds()[i]]);
		}
		vector<double> Gauss_weights;
		vector<Point> Gauss_nodes;
		if (TWO_DIM_LINEAR == basis_type_trial || TWO_DIM_QUADRATIC == basis_type_trial)
		{
			generate_Gauss_local_2D_Tri(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		}
		else if (TWO_DIM_BILINEAR == basis_type_trial || TWO_DIM_BIQUADRATIC == basis_type_trial)
		{
			generate_Gauss_local_2D_Rec(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		}

		//local solution vector
		vector<double> uh_vector;
		for (unsigned int i = 0; i < Tb_test[n].getIds().size(); ++i)
		{
			uh_vector.push_back(solution_vector[Tb_test[n].getIds()[i]]);
		}

		double int_value = Gauss_int_error_2D(analytical_solution, uh_vector, Gauss_weights, Gauss_nodes, vertices, basis_type_trial, der_x, der_y, Tb_test[n].getIds().size());;

		error += int_value;
	}
	return sqrt(error);
}


double compute_L_infi_error(vector<Point> P, vector<Triangle>T, vector<Triangle>Tb_test,
	BASIS_TYPE basis_type_trial, int basis_der_x, int basis_der_y, int Gauss_type, vector<double> solution_vector)
{

	vector<double> r_vec;
	for (unsigned int n = 0; n < Tb_test.size(); ++n)
	{
		vector<Point> vertices;//vertices (coordinates) of nth element(mesh nodes not finite element nodes)		
		for (int i = 0; i < T[n].get_ele_node_num(); ++i)
		{
			vertices.push_back(P[T[n].getIds()[i]]);
		}
		vector<double> Gauss_weights;
		vector<Point> Gauss_nodes;
		if (TWO_DIM_LINEAR == basis_type_trial || TWO_DIM_QUADRATIC == basis_type_trial)
		{
			generate_Gauss_local_2D_Tri(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		}
		else if (TWO_DIM_BILINEAR == basis_type_trial || TWO_DIM_BIQUADRATIC == basis_type_trial)
		{
			generate_Gauss_local_2D_Rec(vertices, Gauss_type, Gauss_nodes, Gauss_weights);
		}
		//local solution vector		
		vector<double> uh_vector;
		for (unsigned int i = 0; i < Tb_test[n].getIds().size(); ++i)
		{
			uh_vector.push_back(solution_vector[Tb_test[n].getIds()[i]]);
		}

		//initial value is the first Gauss_nodes

		double r = analytical_solution(Gauss_nodes[0], 0, 0) - local_FE_function_2D(uh_vector, Tb_test[n].getIds().size(), Gauss_nodes[0], vertices, basis_type_trial, basis_der_x, basis_der_y);
		for (unsigned int k = 0; k < Gauss_nodes.size(); ++k)
		{
			double dr = fabs(analytical_solution(Gauss_nodes[k], 0, 0) - local_FE_function_2D(uh_vector, Tb_test[n].getIds().size(), Gauss_nodes[k], vertices, basis_type_trial, basis_der_x, basis_der_y));
			if (r < dr)
			{
				r = dr;
			}
		}
		r_vec.push_back(r);

	}
	return *max_element(r_vec.begin(), r_vec.end());
}



int main()
{
	double mesh_size = 0.1;
	const double x_min = -1;
	const double x_max = 1;
	const double y_min = -1;
	const double y_max = 1;

	BASIS_TYPE basis_type_trial = TWO_DIM_BILINEAR;//Will be overrided by the element type of the mesh file
	BASIS_TYPE basis_type_test = TWO_DIM_BILINEAR;

	ELE_TYPE ele_type;

	const int basis_der_x_trial = 1;
	const int basis_der_y_trial = 1;

	const int basis_der_x_test = 1;
	const int basis_der_y_test = 1;

	//Gauss quadrature type is 9
	const int Gauss_type = 9;

	cout << "\nHomework of the finite element method course: 2D elliptic equation \n";
	cout << "\nSimulation begins\n";


	int mark;
	string mesh_file_path;
	cout << "\nUse gmsh to generate linear triangle mesh (1) or read mesh from file(2)?\n";
	cin >> mark;
	if (mark == 1)
	{
		cout << "\nInput the mesh size:\n";
		cin >> mesh_size;
	}
	else
	{
		cout << "\nInput the mesh filename\n";
		cin >> mesh_file_path;
	}
	vector<Point> P; vector<Triangle> T;
	vector<Point> Pb_trial; vector<Triangle> Tb_trial;
	vector<Point> Pb_test; vector<Triangle> Tb_test;
	vector<vector<boundary_node>> boundarynodes;//boundary_nodes recorded with each boundary
	vector<vector<boundary_edge>> boundaryedges;//boundary edges recorded with each boundary

	clock_t start, end;
	start = clock();
	if (mark == 1)
	{
		cout << "\ngenerate linear triangle mesh with gmsh" << endl;
		//Generate Pb and Tb with gmsh(Pb=Nodes,Tb=Triangles)
		//also return the bounds nodes and boundary edges
		//this function is only for the linear triangle elements generation
		generate_mesh_with_gmsh(x_min, y_min, x_max, y_max, mesh_size, Pb_test, Tb_test, boundarynodes, boundaryedges, ele_type);

	}
	else
	{
		//read mesh from file generated by gambit
		//also return the boundary nodes and edges
		cout << "read mesh from file" << endl;
		if (!read_mesh_from_file(mesh_file_path, Pb_test, Tb_test, boundarynodes, boundaryedges, ele_type))
		{
			cout << "\n error: read mesh from file failed" << endl;
			return 1;
		}
	}
	end = clock();
	cout << "Time consumption of mesh generation or read:\t" << (double)(end - start) / CLOCKS_PER_SEC << " sec" << endl;

	//output the mesh type
	cout << "\nMesh read or generated success" << endl;
	cout<<" Number of finite element nodes:\t"<<Pb_test.size()<<"\n Number of finite element elements: \t" << Tb_test.size() << endl;
	cout << " Number of each finite element:\t" <<ele_type<< endl;

	//Set the basis function based on the element type which is read from the mesh file
	set_basis_func_type(ele_type, basis_type_trial, basis_type_test);

	//As the mesh generated by gmsh or gambit is the finite element mesh, which means the Pb and Tb not P and T
	//So we need to generate P and T by Pb and Tb
	generate_P_T(basis_type_trial, basis_type_test, P, T, Pb_trial, Tb_trial, Pb_test, Tb_test);

	//assign boundary values
	cout << "\nAssign boundary node and edge values" << endl << endl;
	assign_boundarynodes_value(boundarynodes, Pb_test);

	assign_boundaryedges_value(boundaryedges, Pb_test);


	pair<int, int> matrix_size{ Pb_test.size(),Pb_trial.size() };

	//Assemble matrix
	cout << "Assemble matrix A" << endl;
	start = clock();
	SpMat A = assemble_matrix_2D(matrix_size, P, T, Tb_trial, Tb_test, basis_type_trial, basis_der_x_trial, basis_der_y_trial, basis_type_test, basis_der_x_test, basis_der_y_test, Gauss_type);
	end = clock();
	cout << "Time consumption of matrix assemble:\t" << (double)(end - start) / CLOCKS_PER_SEC << " sec" << endl << endl;

#ifdef _DEBUG
	for (int i = 0; i < A.rows(); ++i)
	{
		for (int j = 0; j < A.cols(); ++j)
		{
			cout << setw(10) << A.coeffRef(i, j) << "\t";
		}
		cout << endl;

	}
#endif

	//Assemble right side vector
	cout << "Assemble right side vector" << endl;
	start = clock();
	Rsv b = assemble_vector_2D(matrix_size, P, T, Tb_trial, basis_type_test, 0, 0, Gauss_type);
	end = clock();
	cout << "Time consumption of right side vector assemble:\t" << (double)(end - start) / CLOCKS_PER_SEC << " sec" << endl << endl;
#ifdef _DEBUG
	for (int i = 0; i < b.size(); ++i)
	{
		cout << b[i] << endl;
	}
#endif

	//treat boundary conditions
	cout << "Treat boundary conditions" << endl;
	start = clock();
	treat_boundary(A, b, boundarynodes, boundaryedges, P, T, Pb_test, Tb_test, basis_type_test,Tb_trial,basis_type_trial, Gauss_type);
	end = clock();
	cout << "Time consumption of boundary conditions treatment :\t" << (double)(end - start) / CLOCKS_PER_SEC << " sec" << endl << endl;

	//Solve the linear equations with Eigen
	//Define a SparseLU solver
	cout << "Solve the linear equations" << endl;
	start = clock();
	Eigen::SparseLU<SpMat> solver;
	solver.compute(A);
	if (solver.info() != Eigen::Success)
	{
		cout << "\nerror: Decomposition failed when using Eigen to compute matrix A\n";
	}
	Rsv x = solver.solve(b);
	if (solver.info() != Eigen::Success)
	{
		cout << "\nerror: Solving failed when using Eigen to solve Ax=b\n";
	}
	end = clock();
	cout << "Time consumption of linear equations solving:\t" << (double)(end - start) / CLOCKS_PER_SEC << " sec" << endl << endl;
	cout << "Simulation end" << endl;

	//Print the solution
	if (x.size() < 500)
	{
		cout << "\nSolution:\n";
		for (int i = 0; i < x.size(); ++i)
		{
			cout << "x[" << i << "]=\t" << setiosflags(ios::scientific) << x(i) << endl;
		}

	}

	//Print the error of the solution at all mesh nodes  
	//Two many output have a low efficiency
	//So if the number of unknown is larger than 500, do not print the solutions
	if (x.size() < 500)
	{
		cout << "\nError of all nodes:\n";
	}
	double max_error = fabs(x(0) - analytical_solution(Pb_trial[0], 0, 0));
	for (int i = 0; i < x.size(); ++i)
	{
		if (x.size() < 500)
		{

			cout << "x[" << i << "]=" << setiosflags(ios::scientific) << fabs(x(i) - analytical_solution(Pb_trial[i], 0, 0)) << endl;
		}
		if (max_error < fabs(x(i) - analytical_solution(Pb_trial[i], 0, 0)))
		{
			max_error = fabs(x(i) - analytical_solution(Pb_trial[i], 0, 0));
		}
	}
	cout << "The maximum absolute error at all nodes:\t" << setiosflags(ios::scientific) << max_error << endl;

	vector <double> solution_vec;
	for (int i = 0; i < x.size(); ++i)
	{
		solution_vec.push_back(x[i]);
	}

	//Compute the L inifite norm erro
	cout << "L∞ error=\t" << compute_L_infi_error(P, T, Tb_test, basis_type_test, 0, 0, Gauss_type, solution_vec) << endl;

	//Compute the L2 norm error and H1 semi-norm error

	cout << "L2 error=\t" << compute_Hs_error(P, T, Tb_trial, basis_type_test, 0, 0, Gauss_type, solution_vec) << endl;
	cout << "H1 error=\t" << compute_Hs_error(P, T, Tb_trial, basis_type_test, 1, 1, Gauss_type, solution_vec) << endl;


	//Output the results to tecplot or vtk
	export_to_tecplot(Pb_trial, Tb_trial, "result.plt", x, basis_type_trial);
	export_to_vtk(Pb_trial, Tb_trial, "result.vtk", x);


	return 0;
}
